Model Results, Bias-corrected survey strata

Fall results

Make a lookup table:

# strata.limits <- as.list(c("AllEPU" = AllEPU, 
#                          "MABGB" = MABGB, 
#                          "MABGBinshore" = MABGBinshore, 
#                          "MABGBoffshore" = MABGBoffshore, 
#                          "bfall" = bfall,
#                          "bfallnot" = bfallnot,
#                          "bfin" = bfin,
#                          "bfinnot" = bfinnot,
#                          "bfoff" = bfoff,
#                          "bfoffnot" = bfoffnot,
#                          "MABGBalbinshore" = MABGBalbinshore,
#                          "MABGBoffshorebigin" = MABGBoffshorebigin))

stratlook <- data.frame(Stratum = c("Stratum_1",
                                      "Stratum_2",
                                      "Stratum_3",
                                      "Stratum_4",
                                      "Stratum_5",
                                      "Stratum_6",
                                      "Stratum_7",
                                      "Stratum_8",
                                      "Stratum_9",
                                      "Stratum_10",
                                      "Stratum_11",
                                      "Stratum_12"),
                           Region  = c("AllEPU", 
                                       "MABGB", 
                                       "MABGBinshore", 
                                       "MABGBoffshore", 
                                       "bfall",
                                       "bfallnot",
                                       "bfin",
                                       "bfinnot",
                                       "bfoff",
                                       "bfoffnot",
                                       "MABGBalbinshore",
                                       "MABGBoffshorebigin"))

Fall Index

Plot individual time series (bias corrected) with standard errors:

splitoutput <- read.csv("pyindex/allagg_fall_500_lennosst_split_biascorrect/Index.csv")

splitoutput <- splitoutput %>%
  left_join(stratlook)

ggplot(splitoutput, aes(x=Time, y=Estimate, colour=Region)) +
  geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_point()+
  geom_line()+
  facet_wrap(~Region, scales = "free_y") +
  guides(colour = guide_legend(ncol=2)) +
  #theme(legend.position = c(1, 0),
   #     legend.justification = c(1, 0))
  theme(legend.position="none")

in2off <- splitoutput %>%
  dplyr::select(Time, Region, Estimate, Std..Error.for.Estimate) %>%
  tidyr::pivot_longer(c(Estimate, Std..Error.for.Estimate), names_to = "Var") %>%
  dplyr::group_by(Var) %>%
  tidyr::pivot_wider(names_from = Region, values_from = value) %>%
  dplyr::mutate(AlbInshore = MABGBalbinshore,
                BlueInshore = bfin,
                BlueOffshore = bfoff,
                OthOffshore = MABGB - (bfoff + bfin + MABGBalbinshore),
                SumMABGB = AlbInshore + BlueInshore + BlueOffshore + OthOffshore) %>%
  dplyr::select(Time, AlbInshore, BlueInshore, BlueOffshore, OthOffshore, SumMABGB, MABGB) %>%
  tidyr::pivot_longer(!c(Time, Var), names_to = "Region", values_to = "value") %>%
  tidyr::pivot_wider(names_from = "Var", values_from = "value")

ggplot(in2off, aes(x=Time, y=Estimate, colour = Region)) +
  geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Fall Prey Index, Mid-Atlantic and Georges Bank")

or as proportions (here proportion of MABGB index).

MABGBprop <- in2off %>%
  #dplyr::filter(Region != "AllEPU") %>%
  dplyr::select(Time, Region, Estimate) %>%
  tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
  dplyr::mutate(AlbInshoreprop = AlbInshore/MABGB,
                BlueInshoreprop = BlueInshore/MABGB,
                BlueOffshoreprop = BlueOffshore/MABGB,
                OthOffshoreprop = OthOffshore/MABGB) %>%
  tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate") %>%
  dplyr::filter(Region %in% c("AlbInshoreprop", "BlueInshoreprop", "BlueOffshoreprop",
                              "OthOffshoreprop"))
  

ggplot(MABGBprop, aes(x=Time, y=Estimate, colour = Region)) +
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Fall Prey Index as proportion of Mid-Atlantic and Georges Bank")

Fall predicted ln-density

Fall density maps with covariates

Fall Diagnostics

diagplots <- c("Data_and_knots",
               "Data_by_year",
               "quantile_residuals",
               "quantile_residuals_on_map",
               "Aniso",
               "center_of_gravity")

for(p in diagplots){
  
    cat("  \n####",  as.character(p),"  \n")
    cat(paste0("![](pyindex/allagg_fall_500_lennosst_split_biascorrect/",
                      p,
                      ".png)")) 
    cat("  \n")   
    
  }

Data_and_knots

Data_by_year

quantile_residuals

quantile_residuals_on_map

Aniso

center_of_gravity

Spring results

Spring Index

Plot individual time series

splitoutput <- read.csv("pyindex/allagg_spring_500_lennosst_split_biascorrect/Index.csv")

splitoutput <- splitoutput %>%
  left_join(stratlook)

ggplot(splitoutput, aes(x=Time, y=Estimate, colour=Region)) +
  geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_point()+
  geom_line()+
  facet_wrap(~Region, scales = "free_y") +
  guides(colour = guide_legend(ncol=2)) +
  #theme(legend.position = c(1, 0),
   #     legend.justification = c(1, 0))
  theme(legend.position="none")

or just the indices from inshore (alb), inshore bluefish, offshore bluefish, and further out

in2off <- splitoutput %>%
  dplyr::select(Time, Region, Estimate, Std..Error.for.Estimate) %>%
  tidyr::pivot_longer(c(Estimate, Std..Error.for.Estimate), names_to = "Var") %>%
  dplyr::group_by(Var) %>%
  tidyr::pivot_wider(names_from = Region, values_from = value) %>%
  dplyr::mutate(AlbInshore = MABGBalbinshore,
                BlueInshore = bfin,
                BlueOffshore = bfoff,
                OthOffshore = MABGB - (bfoff + bfin + MABGBalbinshore),
                SumMABGB = AlbInshore + BlueInshore + BlueOffshore + OthOffshore) %>%
  dplyr::select(Time, AlbInshore, BlueInshore, BlueOffshore, OthOffshore, SumMABGB, MABGB) %>%
  tidyr::pivot_longer(!c(Time, Var), names_to = "Region", values_to = "value") %>%
  tidyr::pivot_wider(names_from = "Var", values_from = "value")

ggplot(in2off, aes(x=Time, y=Estimate, colour = Region)) +
  geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Spring Prey Index, Mid-Atlantic and Georges Bank")

or as proportions (here proportion of MABGB index).

MABGBprop <- in2off %>%
  #dplyr::filter(Region != "AllEPU") %>%
  dplyr::select(Time, Region, Estimate) %>%
  tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
  dplyr::mutate(AlbInshoreprop = AlbInshore/MABGB,
                BlueInshoreprop = BlueInshore/MABGB,
                BlueOffshoreprop = BlueOffshore/MABGB,
                OthOffshoreprop = OthOffshore/MABGB) %>%
  tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate") %>%
  dplyr::filter(Region %in% c("AlbInshoreprop", "BlueInshoreprop", "BlueOffshoreprop",
                              "OthOffshoreprop"))
  

ggplot(MABGBprop, aes(x=Time, y=Estimate, colour = Region)) +
  geom_point()+
  geom_line()+
  #facet_wrap(~Region) + #+ , scales = "free_y"
  #theme(legend.position = c(1, 0),
  #      legend.justification = c(1, 0))
  ggtitle("Spring Prey Index as proportion of Mid-Atlantic and Georges Bank")

Spring predicted ln-density

Spring density maps with covariates

Spring Diagnostics

diagplots <- c("Data_and_knots",
               "Data_by_year",
               "quantile_residuals",
               "quantile_residuals_on_map",
               "Aniso",
               "center_of_gravity")

for(p in diagplots){
  
    cat("  \n####",  as.character(p),"  \n")
    cat(paste0("![](pyindex/allagg_spring_500_lennosst_split_biascorrect/",
                      p,
                      ".png)")) 
    cat("  \n")   
    
  }

Data_and_knots

Data_by_year

quantile_residuals

quantile_residuals_on_map

Aniso

center_of_gravity

All results in the respective pyindex/allagg_fall_500_lennosst_split_biascorrect and allagg_spring_500_lennsst_splitbiascorrect folders.

The full results are on google drive rather than github to save space.

Still to do:

  • fix index within 3 miles of shore and outside that (highest priority)

References

Garrison, L., and Link, J.
2000. Dietary guild structure of the fish community in the Northeast United States continental shelf ecosystem. Marine Ecology Progress Series 202:231–240. https://doi.org/10.3354/meps202231.
Ng, E.L., Deroba, J.J., Essington, T.E., Grüss, A., Smith, B.E., and Thorson, J.T.
2021. Predator stomach contents can provide accurate indices of prey biomass. ICES Journal of Marine Science 78:1146–1159. https://doi.org/10.1093/icesjms/fsab026.
Perretti, C.T., and Thorson, J.T.
2019. Spatio-temporal dynamics of summer flounder (Paralichthys dentatus) on the Northeast US shelf. Fisheries Research 215:62–68. https://doi.org/10.1016/j.fishres.2019.03.006.
Schoener, T.W.
1970. Nonsynchronous Spatial Overlap of Lizards in Patchy Habitats. Ecology 51:408–418. https://doi.org/10.2307/1935376.
Smith, B.E., and Link, J.S.
2010. The Trophic Dynamics of 50 Finfish and 2 Squid Species on the Northeast US Continental Shelf. NOAA Technichal Memorandum NMFS-NE-216. National Marine Fisheries Service, 166 Water Street, Woods Hole, MA 02543-1026.
Thorson, J.T.
2019. Guidance for decisions using the Vector Autoregressive Spatio-Temporal (VAST) package in stock, ecosystem, habitat and climate assessments. Fisheries Research 210:143–161. https://doi.org/10.1016/j.fishres.2018.10.013.
Thorson, J.T., and Barnett, L.A.K.
2017. Comparing estimates of abundance trends and distribution shifts using single- and multispecies models of fishes and biogenic habitat. ICES Journal of Marine Science 74:1311–1321. https://doi.org/10.1093/icesjms/fsw193.
Thorson, J.T., and Kristensen, K.
2016. Implementing a generic method for bias correction in statistical models using random effects, with spatial and population dynamics examples. Fisheries Research 175:66–74. https://doi.org/10.1016/j.fishres.2015.11.016.